1. Field of the Invention
This invention pertains to the data processing of seismic traces to suppress noise and more specifically to such processing of stacked data through the development and employment of improved weighting factors.
2. Description of the Prior Art
This invention is relevant to acquisition of seismic data in the field of geophysical exploration for petroleum and minerals. Although it is not restricted to exploration in either land or marine environments, this method is presented with primary emphasis on land seismic applications since this is the area in which it has the most direct and widespread application.
The general seismic prospecting method involves transmission of acoustic waves through the earth and reception of these waves at the earth's surface via transducers called geophones. Typically, these geophones are wired in parallel collections called groups in order to increase electrical signal levels and attenuate surface traveling noises. The acoustic waves may be generated by various types of sources, dynamite and hydraulic vibrators being the most common. As these waves propagate downward through the earth, portions of their energy are sent back toward the earth's surface through the effects of reflections and refractions which occur whenever abrupt changes in acoustic impedance are encountered. Since these impedance changes coincide with sedimentary layer boundaries it is possible to image the layers by appropriate processing of the signals returned to geophone groups.
The quality of the images obtained from seismic data processing is limited by three principal factors. First, the factor of spherical divergence causes the signal amplitude to decay in a manner which is inversely proportional to the distance traveled. Second, the effects of absorption due to the inelasticity of the media of propagation further reduces the available signal level. Finally, the effects of source generated and environmentally generated noise tend to mask whatever signal remains after divergence and absorption effects. In fact, it is the effect of noise which limits the depth to which the seismic method can penetrate. To some extent this limitation can be overcome by increasing the source signal strength so that the received signal level is raised relative to the noise. However, since there are practical limitations on the size of the source, a preferable approach is to suppress the noise level as much as possible.
The amount of achievable noise reduction and the techniques used to suppress noise differ according to the type of noise being attacked. Since source generated noise is not random and is repeatable over various trials it is difficult if not impossible to suppress it until the data is processed. On the other hand, noise which is not source generated is generally random and, even if it is not, it is almost certainly asynchronous with the data gathering sequence. These factors allow environmentally generated noise to be suppressed by techniques performed as the data is being acquired.
There are three techniques which are effective at reducing environmentally generated noise. First, it is possible in some cases to detect and edit noise bursts out of the acquired data. This technique is therefore referred to as noise or trace editing. Second, it is possible to cancel random noises by summing like traces acquired during successive repetitions of the data acquisition sequence. This process is referred to as vertical stacking. Finally, the process of correlation provides an improvement in signal-to-noise ratio in cases where vibratory sources are used.
Seismic sources fall within one of two classifications, impulsive and non-impulsive. Impulsive sources are those which produce energy distributions which are highly time-limited but whose energy concentration is high. Typically, the signals produced by these are considered minimum phase or near-minimum phase time series, indicating that the rate of energy generation is as fast as possible for a given spectral energy distribution. As the bandwidth of such sources increases their corresponding output signals, or signatures, collapse in time. In the limit of infinite bandwidth, such signatures would reduce to what is referred to in signal processing literature as a delta function, a pulse of infinitesimal width and possessing infinite amplitude. Such pulses typically require little or no processing to reduce the pulse duration to acceptable limits. Conversely, non-impulsive sources are characterized by waveforms which are distributed in time which can have spectral energy densities which are comparable to impulsive sources but whose signatures persist in time and with lower amplitude. The signals produced by such sources are termed mixed-phase since the rate of energy generation is not as fast as their corresponding spectra will allow. One common conception of signals from non-impulsive sources is to consider them as concatenated collections of minute impulses of varying magnitude and sign.
Data collected from the use of non-impulsive sources such as vibrators always requires processing aimed at reducing the effective time duration of the signature. In such cases, a process known as correlation or matched filtering is applied to the data to collapse the pulse primarily by rotating the phase spectrum of the associated signature. While this process requires additional steps in processing, it does provide a certain amount of noise rejection. It also allows certain trace editing procedures to be performed prior to correlation which serve to suppress noise energy without seriously impacting the signal content of the processed output. This is intuitively evident by observing first of all that since a non-impulsive signal is broadly distributed in time, it is possible to remove small segments of it without detracting significantly from the total signal energy. Second, it can be observed that correlation applies the same phase adjustments to noise energy such that localized noise bursts which were once confined to a small time segment are spread over a broad segment of the data after correlation. This later principle has the effect of increasing the signal-to-noise ratio of the correlated data provided that the noise bursts are not large and are randomly distributed throughout the traces.
In considering the effect of noise upon the output of the correlation procedure, it is first necessary to define two classifications of noise. The first type of noise to be considered is that which persists in time with the same statistical characteristics, such as noise associated with certain types of stationary mechanical and electrical equipment as well as ground unrest due to other noise generators that are asynchronous in nature and are far removed so that no single noise source is dominant. In statistical terms, such noise sources are characterized by probability density functions which are constant or "stationary" with time and as such they are termed stationary noise sources. Conversely, the second classification of noise sources involves those that are sporadic in nature, such as lightning and activity of the crew along a seismic line. Such noise sources are characterized by probability density functions that are time-variant and are termed "non-stationary" noise sources.
Differentiation between stationary and non-stationary noise sources can be ambiguous since the probability density functions obtained vary considerably with the size of the time gate over which the noise source is analyzed. For example, a noise source which produces pulses at half-second intervals would have a constant probability density function and would appear to be stationary if the analysis time gate were sufficiently long. Conversely, if the time gate is small the estimate of the probability density function tends to vary according to the number of noise pulses each contains and in such cases the noise would appear to be non-stationary. In the seismic case, the issue of stationarity is generally resolved by using time gates which are comparable to the period of the lowest frequency included in the signal.
The relative effects of stationary and non-stationary noises on correlated data quality can be seen by considering the nature of these noises and their effect on the correlation process. The phase shifts that occur during the correlation process have no effect on stationary noises since by definition stationary noises are already distributed throughout each input trace. Non-stationary noise sources on the other hand, can and typically do produce noise energy that arrives in bursts. The phase shifts associated with correlation then effectively smear the energy contained in the bursts into surrounding time segments. This will have effects that vary from being unnoticeable to that of corrupting the correlated trace and possibly the output trace into which it is summed.
The effect of correlation on data that contains noise bursts can be handled more analytically by recognizing that correlation is quite similar to convolution, the only difference being that in correlation one of the operand traces is effectively time reversed. By analogy, a noise burst embedded in an input trace is similar to a filter impulse response which is convolved with an input time series. The net effect on the output is to overlay a filtered and time reversed version of the correlation reference trace on top of the correlated trace. By this observation it is apparent that not only the amplitude of the noise burst but also its shape have a bearing on the manner in which the noise is transmitted into the correlated trace.
Considerations concerning the shape of noise bursts point up some potential pitfalls of noise limiting techniques since they can have a significant effect on the shape of noise bursts. Consider, for example, the pulse which results when an impulse is applied to a inelastic medium. This is very appropriate since the surface of the earth is generally very inelastic. The initial compressional lobe of the pulse peaks at a large amplitude are followed by a relaxational lobe that has substantially the same area but a lower peak amplitude and longer duration. The approach taken by some traditional limiter techniques would be to detect and clip the initial lobe but leave the second lobe intact since it is of lower amplitude. This approach has the intended effect of reducing the peak-to-peak amplitude of the pulse, however it also modifies the spectrum of the noise burst by introducing an additional low frequency content that could still correlate well with portions of the reference signal. In such cases, the net effect of clipping may be trade off the frequencies at which the noise burst correlates with the reference. In light of this, the ideal solution to the limiter problem would be to somehow detect and scale the entire noise burst down in amplitude to some optimal value.
One problem associated with the elimination of noise bursts is first, that it is difficult to determine relative amounts of signal and noise, and second, that it is difficult to discriminate against noise without also affecting signal. In the example given above, the only reason that the noise burst was recognized as such was that it has an anomolously high amplitude. On a trace by trace basis without any information other than the expected level and approximate bandwidth of the signal, the only criteria for noise burst detection is the presence of obvious amplitude anomalies. One advantage of using non-impulsive sources is that due to their mixed-phase quality their signal strength is relatively low and sharp changes in signal strength are rare. It is, therefore, easier to detect noise bursts due to their contrasting qualities of high levels and sharp onset.
Assuming that a noise burst is first detected, the problem then arises concerning the effective elimination of the noise while minimizing any effects on signal. Furthermore, as stated above, the problem of eliminating the noise burst without causing problems in the correlation results also exists. In addition to the basic technique of clipping noise bursts, several more elaborate schemes have been developed. A typical example of this is known as the recursive edit in which a time segment of the trace is zeroed starting at the last zero crossing before the noise burst and lasting for a specified amount of time after the burst. By zeroing all data between zero crossings this method attempts to minimize generation of high frequencies as well as the edge effects on correlation. Clearly, however, this does not recover the signal which was in the time segment along with the noise burst.
A method first employed by Embree and discussed in his U.S. Pat. No. 3,398,396, neatly addresses problems of noise detection and suppression in one straightforward procedure. This procedure is called the Diversity Stack and has found wide application in suppressing non-stationary noise while preserving signal.
The power of the Diversity Stack procedure comes from exploitation of the fact that signal is invariant between input traces originating from the same geodetic receiver location and relationship to the source location. Assuming that signal is uncorrelated with noise the total energy within a trace time segment is given by: ##EQU1## (where t1, t2=begin and end times of time segment.)
If it is assumed that S(t) is invariant between like trace time segments, any relative increase in Pi must be due to an increased noise level for that trace segment. Using this fact, the Diversity Power Stack prorates corresponding segments of input traces according to the inverse of the energy within that segment before stacking each trace. A weight stack is also maintained as a record of the total weight applied to each trace sample so that the aggregate weighting can be removed when the stack is complete. This last step is referred to as stack normalization.
The following is a formal detailed analysis of optimal stacking in the presence of time-varying noise components which are characterized by sporadic increases in noise amplitudes:
Given a set of trace time series which originate from a single surface location as the result of multiple seismic excitations at a substantially constant surface location, each trace within the set may be described by G.sub.i (t)=S.sub.i (t)+N.sub.i (t) where S.sub.i (t) is the signal component containing useful seismic information and N.sub.i (t) is the noise component consisting of sporadic time-variant noises and ambient time-invariant noise. The subscript .sub.i is used to indicate the trace within the set which arises from the i'th seismic excitation at the source location. A generalized weighted trace stack is then given by ##EQU2##
By definition the signal S.sub.i (t) is independent of i so that the above may be rewritten as ##EQU3##
The total energy in G(t) is ##EQU4##
Here the first term is the total signal energy and the second term represents the total noise energy in G(t) after weighting, stacking and normalization.
The optimal value for an arbitrary weight W.sub.k is obtained by minimizing the noise energy term with respect to W.sub.k. Differentiating with respect to W.sub.k ##EQU5##
Setting the numerator to zero and simplifying the following relationship is obtained: ##EQU6##
Upon inspection, it is apparent that the relationship W.sub.i =C/N.sub.i for i=1,M provides a general solution where C is an arbitrary constant.
Although W.sub.i is not time varying, the above derivation may be extended to show that W.sub.i (t)=1/N.sub.i (t) is an optimal time varying weight function where N.sub.i (t) is the noise energy in a time gate containing t. The only restriction is that this time gate be large enough to obtain an accurate estimate of N.sub.i (t).
Substituting W.sub.i =W.sub.i (t), S=S(t) and N.sub.i =N.sub.i (t) into the above equations EQU W.sub.i (t)=1/N.sub.i (t) for optimal W.sub.i (t) ##EQU7##
Assuming optimal weighting N(t) becomes ##EQU8##
The signal to noise power ratio is then ##EQU9## where .theta..sub.i (t) is the signal to noise power ratio of the component traces. The signal to noise amplitude ratio is ##EQU10##
From these relationships for .theta.(t) and .phi.(t), it is evident that the signal to noise power ratio of the optimally weighted stack is equal to the sum of the signal to noise power ratios of the component traces. Furthermore when N.sub.i (t), i=1, M are all stationary random time series with equal variance the signal to noise amplitude ratio reduces to EQU .phi.(t)=M.sup.1/2 S(t).sup.1/2 /N(t).sup.1/2
Thus, the improvement in signal to noise ratio of the weighted stack is equivalent to that obtainable by unweighted stacking. In fact since W.sub.i (t)=1/N.sub.i (t)=1/N, the weighted stack becomes an unweighted stack in the absence of non-stationary noise components.
Since the actual relationship of signal and noise energies is not generally known, the approximation for W.sub.i (t) must be based upon (S(t)+N.sub.i (t)) instead of N.sub.i (t) alone. The first and most obvious approximation takes the form W.sub.i '(t)=C'/(.delta.(t)+N.sub.i (t)), but in this case W.sub.i '(t) underestimates W.sub.i (t) except when S(t)&lt;&lt;N.sub.i (t).
In order to overcome this underestimation and optimize W.sub.i (t) for estimated signal to background noise ratios an additional parameter B.sub.i (t) is added so that the approximation becomes EQU W'(t)=C'/(.phi.(t)+N.sub.i (t)).sup.Bi(t)/2
where B.sub.i (t).gtoreq.2
Since only relative weighting is of concern, the objective in evaluating B.sub.i (t) is not to force W.sub.i '(t)=W.sub.i (t) but rather dW.sub.i '(t)/dN.sub.i (t)=dW.sub.i (t)/dN.sub.i (t), or alternatively EQU d Ln(W.sub.i '(t))/dN.sub.i (t)=d Ln(W.sub.i (t))/dN.sub.i (t)
Differentiating and solving for B.sub.i (t) gives EQU d Ln(W.sub.i '(t))/dN.sub.i (t)=-B.sub.i (t)/2(N.sub.i (t)+S(t)) EQU d Ln(W.sub.i (t))/dN.sub.i (t)=-1/N.sub.i (t) EQU 1/N.sub.i (t)=B.sub.i (t)/2(N.sub.i (t)+S(t)) EQU B.sub.i (t)=2(1+(S(t)/N.sub.i (t))
This relationship allows B.sub.i (t) to be calculated based upon estimated or anticipated signal to noise ratios. Theoretically, there should be no upper bound on B.sub.i (t). From a practical standpoint, however, it can be argued that when S(t)/N.sub.i (t)&gt;&gt;1, the signal to noise ratio of the component traces G.sub.i (t) is already so large that any further increase in B.sub.i (t) and consequently the weight function W.sub.i (t) would have only marginal value. It may also be observed that when S(t)&lt;&lt;N(t), B(t) reduces to the value 2.0 as embodied in the Diversity Power Stack, consistent with the fact that it also assumes S(t)&lt;&lt;N(t).
As mentioned previously the signal to noise power ratio S(t)/N.sub.i (t) is not generally known unless some type of coherence analysis is performed prior to stacking, but to do this on a large scale would certainly be cost prohibitive. Since the purpose of B.sub.i (t) is to equivalence noise dependent variations in the optimal weight W.sub.i (t) and the approximate optimal weight W.sub.i '(t), some efficient criteria for estimating B.sub.i (t) is needed in lieu of coherence analysis procedures. Here, a fairly accurate estimate can be made of expected noise power using time gates at the start and end of one or more component traces and assuming that the power observed in these gates is solely due to noise with both stationary and non-stationary components which are assumed to be uncorrelated. This estimate is given by EQU Ne=median (N.sub.ij), i-1,M; j=1,N
where N.sub.ij is the noise estimate associated with the j'th noise estimation gate of the i'th input trace.
The justification in using the median as opposed to another estimator such as the average is that it is not overly influenced by large noise burst which happen to fall within an estimation gate. This tends to give an estimate for Ne that is nested within the noise estimates of the majority of the time gates.
Using this estimate the following approximations are made for S(t) and B.sub.i (t): ##EQU11##
This last expression allows the function B.sub.i (t) to be computed on a per trace basis using quantities that are directly obtainable from the data, the total trace energy G.sub.i (t) and the background noise estimate N.sub.e.
The lower bound on B.sub.i (t) is set due to knowledge that optimal B.sub.i (t) never goes below 2 and values of B.sub.i (t) which violate this constraint must be due to overestimation of Ne. By setting the lower limit of B.sub.i (t) at 2, the risks associated with overestimation of the noise power are minimized by restricting B.sub.i (t) from assuming values which are known to be sub-optimal in any case.
Since Embree's development of the Diversity Power Stack, Robinson, "Statistically Optimal Stacking of Seismic Data", Geophysics, Jun. 1970, and Smith et al., U.S. Pat. No. 4,561,075, have also addressed the theory of optimal stacking, both of these treatments further refined the development of the weight function. Methods proposed and evaluated by Robinson were directed toward optimizing weight functions with a minimal number of assumptions about the relative levels of signal and noise energies within the constituent traces. For this reason, the methods of weighting that Robinson developed were very robust but required a large amount of computation and memory storage capacity which made them impractical for real time field applications. Smith et al. on the other hand set forth a recording system which used the basic prior art Diversity Stacking method but went further in allowing selection of integer exponents other than 2 to compensate for signal to noise ratios that are comparable to or in excess of 1. The selection of the exponent was then based upon the expected relationship of signal and noise typical among all traces. Exponents which are greater than two result in a more dramatic suppression of trace amplitude for a given increase in noise energy than do exponents of lower orders. This provides a more optimal weighting of traces during periods of high signal to noise ratio with the consequence that later portions of the traces where the signal component has decayed are overly suppressed in the presence of sporadic noise bursts.
Traditional methods of computing the Diversity Stack typically use time domain linear interpolation between time segments to approximate the weight function, W(t), which causes two fundamental problems. First, the derivatives of W(t) are discontinuous, thus producing high frequencies in the spectrum of W(t). Since multiplication in the time domain is equivalent to convolution in the frequency domain, these higher frequency components will cause spectral smearing in the weighted trace. The second problem is that these methods produce inaccurate estimates of W(t), since the time relationship between the noise bursts and time segments is not controllable.
Therefore, it is a feature of the present invention to provide the development of an improved optimized seismic stack trace utilizing a weighted trace stack and a weight stack wherein the individual weight functions are determined on a one-to-one correspondence with the individual seismic traces following the removal of gross amplitude anomalies, the resulting weighted trace stack then being divided by the resulting weight stack.
It is another feature of the present invention to provide an improved procedure for developing the weights for individual traces whereby the development employs conversions of individual sampled digital data from the time domain to the frequency domain and back again to minimize the mathematical manipulations required for application of filters and the time decimation and converse resampling of time domain counterparts.